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classical Boltzmann equation in the quantum collision regime 
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In this paper we develop a direct simulation Monte Carlo (DSMC) method for simulating highly nonequi- 
librium dynamics of nondegenerate ultra cold gases. We show that our method can simulate the high-energy 
collision of two thermal clouds in the regime observed in experiments [Thomas et al. Phys. Rev. Lett. 93, 1 7320 1 
(2004)], which requires the inclusion of beyond s-wave scattering. We also consider the long-time dynamics 
of this system, demonstrating that this would be a practical experimental scenario for testing the Boltzmann 
equation and studying rethermalization. 
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I. INTRODUCTION 

Within ultra-cold-atom research, there are a range of prob- 
lems requiring the understanding of the dynamics of a normal 
gas. For example, studies of collective modes of Bose ^ and 
Fermi [fS] gases (also see Refs. ||3l01), spin waves Sla], hy- 
drodynamic expansion of a Bose gas near the critical temper- 
ature 101, and more recently, the dynamics and thermalization 
of a nearly degenerate gas of polar molecules HI. These are 
all regimes in which the Boltzmann equation is thought to pro- 
vide an accurate description. In many of these cases, the sys- 
tem is only weakly disturbed from equilibrium, and some ap- 
proximate solution can be provided using a relaxation approx- 
imation for the collision integral and some form of lineariza- 
tion 1^, scaling iflolfTTll . or variational II 1211 ansatz. For more 
strongly dynamical situations, these approaches are insuffi- 
cient, however, the direct solution of the Boltzmann equation 
for the six-dimensional distribution function is generally con- 
sidered intractable and is normally tackled using some form of 
stochastic particle simulation. Some applications of such cal- 
culations include the work of Wu and co-workers lUsUTsll on 
evaporative cooUng and expansion dynamics, Jackson and co- 
workers |[l6l - l20ll on bosonic collective-mode dynamics (cou- 
pled to a superfluid by the Zaremba-Nikuni-Griffin (ZNG) 
formalism uM ). the work of Urban and Schuck j22|. Urban 
ll23i I24I1 . and Lepers et al. Izsll in formulating fermion dy- 
namics (see also Refs. ll26l429ll ). and Barletta et al. ifsoll and 
Barletta llsill in describing sympathetically cooled molecular 
gases. 

Here, we develop an algorithm for simulating the Boltz- 
mann equation that is significantly more accurate and efficient 
than these previous methods and is applicable to more extreme 
regimes of dynamics. Indeed, our main motivation was to de- 
velop a theory capable of describing the ultra-cold-atom col- 
lider developed by the Otago group ||32| - |34|| . In those experi- 
ments (nonquantum degenerate), clouds of bosonic atoms at a 
temperature of ^ 200 nK were accelerated and were collided 
at an energy of ^ 200 /iK (see Fig. [1]). Several features of 
these experiments make the numerical simulation difficult: 

(i) The system is far from equilibrium and accesses a large 
volume of phase space. A good representation of each 
cloud before the collision requires nano-Kelvin energy 
resolution, however, during the collision, atoms are 
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FIG. 1 . (Color online) Ultra cold atom collider: (a) schematic of the 
precollision arrangement of two clouds at ~ 200 nK approaching 
at a collision energy of ^ 200 /iK; (b) schematic of a postcollision 
system, (c) and (d) experimental images of post scattering density 
for two collision energies spanning the d-wave shape resonance, (e) 
and (f) show the theoretical calculations matching the experimental 
results using the direct simulation Monte Carlo (DSMC) method de- 
veloped in this paper. 



scattered over states on the collision sphere with an en- 
ergy spread on the order of a milli-Kelvin. 

(ii) The collision energies are sufficiently large that an ap- 
preciable amount of higher-order (i.e. beyond s-wave) 
scattering occurs. In particular, in experiments p-wave 
scattering ll34ll and a d-wave ll32ll shape resonance have 
been explored. 
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The algorithm we develop is suitable for this regime, and, as 
shown in Figs. [nc)-[nf), it can provide a quantitative model 
for the experimental data in Ref. [l32ll . Feature (i) discussed 
above presents a great challenge, and using the traditional 
Boltzmann techniques employed to date in ultra-cold-atom re- 
search, this would require super computer resources. We show 
how to make use of an adaptive algorithm (that adapts both 
the spatial grid and the times steps to place resources where 
needed) to accurately simulate an ultra cold collider on com- 
modity personal computer hardware. 

We note that, in addition to collider experiments, a capable 
Boltzmann solver would allow theoretical studies in a range of 
areas of emerging interest, such as the turbulence and flow in- 
stabilities in the normal phase of a quantum gas. Here, we will 
focus on the classical regime where the phase-space density is 
small compared to unity such that the many-body effects of 
Bose-stimulated or Pauli-blocked scatterings are negligible. 
However, the systems we consider will be in the quantum col- 
lision regime, whereby the thermal de Broglie wavelength is 
larger than the typical range of the interatomic potential. No- 
tably, in this regime, the scattering is wave like, and quantum 
statistics on the two-body level gives rise to profound effects 
in the individual collision processes, even though many-body 
quantum statistics is unimportant. 

All of the Boltzmann simulations appearing in the ultra- 
cold-atom literature have been based on DSMC-like methods, 
typically employing the algorithm described in Bird's 1994 
monograph IBSll . However, a challenging feature of ultra cold 
gases is that the local properties (e.g., the density) can vary by 
orders of magnitude across the system, and no single global 
choice of parameters for the DSMC can provide a good de- 
scription across this entire range. For this reason, we intro- 
duce the use of two locally adaptive schemes to allow the sys- 
tem to refine the description and to allocate more computa- 
tional resources to regions of high density. These schemes, 
which we discuss in Sec. |III1 are as follows: locally adaptive 
time steps (LATSs) and locally adaptive cells (LACs). 

In Sec. lIVI we validate our algorithm using a variety of tests 
to demonstrate its applicability and performance. Then, in 
Sec. |V] we apply it to the regime of the ultra-cold-atom col- 
lider experiments 1I32I1 . 



11. THEORY 

A. Boltzmann equation 

The system is described semiclassically by the phase space 
distribution function / = / (r, r, t), which evolves according 
to the Boltzmann equation Qql . 



d 
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VrU{r,t) ■ Vp 



/ = /[/], (1) 



and the position space density of the atoms n (r, t) is given by 



The left-hand side of Eq. ([TJ describes the evolution of atoms 
under the potential U (r, t). In general, U (r, t) may contain 
a mean-field term, however, for our analysis in this paper, we 
only consider the case where U (r, t) is an external trapping 
potential. 

The collision integral / [/], accounts for the collisions be- 
tween atoms and is given by 



I[.f] = 



d^Pi 

/l3 



dn^ |pi - p| [/'/{ - //i] , (3) 



where ^ is the differential cross section and /i = 
/ (pi, r, t), f = f (p', r, t), etc. When considering the flow 
of atoms through phase space due to collisions, / [/] has a 
simple interpretation. The term in Eq. (|3) containing ffi de- 
scribes collision events where the atoms are initially at the 
phase-space points (p,r) and (pi,r) and have final states 
(p', r) and {p[, r). The rate of such a colUsion depends on 
the densities of the initial states //i, kinetic factors described 
by the differential cross section, and the flux of incident parti- 
cles, which is proportional to |pi — p|. The opposite process 
where atoms scatter to (p, r) and (pi, r) is accounted for by 
the /'/{ term. The quantum statistics of the atoms can be 
included by the addition of (1 ± /') (1 ± /() terms in the col- 
lision integral, which account for Bose-stimulated scattering 
(+) or Pauli blocking (— ). Here, we neglect quantum statis- 
tics as is appropriate for the regime where / <C 1 and will 
address considerations for the full quantum Boltzmann equa- 
tion elsewhere OTIl . 



B. Partial-wave treatment of collisions 

While our interest here is in ultra cold gases with suffi- 
ciently low phase-space density to neglect many-body quan- 
tum statistics, the two-body collisions themselves are in the 
quantum collision regime and are conveniently characterized 
in terms of a partial- wave expansion. The differential cross 
section for identical bosons (+) or fermions (— ) in the same 
internal state is 



where 



fsc (O) 



h ^ 



0)r 



imv. 



- ^ (2/ + 1) {e^'^' 



cos ( 



(4) 



(5) 



(2) 



is the scattering function, Vr is the magnitude of the relative 
velocity of the colliding particles, 61 is the phase shift associ- 
ated with partial wave /, Pi (cos 9) is the Legendre poly- 
nomial, and 6 is the centre-of-mass scattering angle. In gen- 
eral, the phase shifts have a collision energy dependence (vr), 
which is a nontrivial task to calculate. 

For bosons (fermions), the total wave function is required to 
be symmetric (antisymmetric), and hence, only the even (odd) 
I terms in Eq. ^ contribute to the differential cross section. 
In this paper, we focus on the case of bosons, motivated by the 
experimental work we seek to describe ll32ll . 
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III. DSMC METHOD 
A. Background for DSMC 

The DSMC method is the most widely used tool for model- 
ing fluid flow on the subcontinuum scale and has found itself 
successfully applied to a huge range of physics from shock 
waves ifssll and Rayleigh-Benard flow OSll to aerodynamics 
of spacecraft llsoll . chemical reaction s llioll . microfluidics iiUl . 
acoustics on Earth, Mars, and Titan ll42ll . volcanic plumes on 
Jupiter's moon lo ll43ll . and much more. 

These situations are characterized by being dilute (two- 
body collisions) and having a high Knudsen number (Kn), 
which is given by the ratio of the mean-free path to the rep- 
resentative length scale of the system. For Kn > 0.1, a mi- 
croscopic kinetic theory is necessary, while for Kn < 0.1, 
the system tends to be sufficiently hydrodynamic for a con- 
tinuum approach to be applicable for understanding coarse- 
grained dynamics^ This is not to say that DSMC is inappli- 
cable or is inefficient in this regime; indeed, recently. Bird 
has shown that, in nonequilibrium situations with Kn ~ 0.01, 
the DSMC algorithm (employing many of the techniques we 
introduce for cold atoms here) can be more accurate and ef- 
ficient than Navier-Stokes methods, while also providing de- 
tails of the microscopic (subcontinuum) dynamics MM . We 
also note that the consistent Boltzmann algorithm 14511 was 
developed by making an adjustment to the DSMC algorithm, 
where the positional shifts associated with collisions are taken 
into account, giving the correct hard-sphere virial. This allows 
for exploration into even lower Kn and has been explored in 
the context of quantum nuclear flows ll46ll47ll . 

For reference, cold-atom experiments often operate in the 
collisionless regime (Kn > 1), however, values of Kn ~ 0.01 
have been explored, e.g., the above-critical temperature col- 
lective modes of a ^■^Na gas studied by Stamper- Kum et al. 10] 
had Kn ^ 0.1; Shvarchuck et al. iQl] studied the hydrodynam- 
ical behavior of a normal *^Rb gas in which Kn ~ 0.02 — 0.5. 

B. Overview of formalism and general considerations 

In the DSMC method, the distribution function is repre- 
sented by a swarm of test particles, 

f (p, r,t)^ah^J2^^P- (^)] t'" - (^)] ' 

i=l 

where a = A/p /Mt is the ratio of physical atoms (Afp) to test 
particles (N't)- These test particles are evolved in time in such 
a manner that / (p, r, t) evolves according to the Boltzmann 
equation. 

The basic assumption of DSMC is that the motion of atoms 
can be decoupled from collisions on time scales much smaller 



In the cold-atom community, it is more common to specify tliese regimes 
as uiT, where oj is the excitation frequency and r is the collision time. 



than the mean-collision time. In practice, this means that 
a simulation is split up into discrete time steps At, during 
which, the test particles undergo a collisionless evolution, then 
collisions between test particles are calculated. 

The relation of the test particles to physical atoms is appar- 
ent in Eq. (|6]l when a = 1, but, in general, they are simply a 
computational device for solving the Boltzmann equation. In 
many conventional applications of DSMC, good accuracy can 
be obtained with a S> 1 (i.e., each is a super particle repre- 
senting a larger number of physical atoms), however, in our 
applications on nonequilibrium dynamics of ultra cold gases, 
we often require a <C 1. Increasing the number of particles 
improves both the accuracy and the statistics of the simula- 
tion, and in highly nonequilibrium situations, it can be essen- 
tial to have a large number of particles. The DSMC method 
is designed so that the number of computational operations 
per time step scales linearly with the number of particles, i.e., 
O (A/t)- The recent work of Lepers et al. Izsll departs from 
DSMC by using a stochastic particle method similar to that 
developed in nuclear physics for the simulation of heavy-ion 
collisions ll48[ l49ll . which tests if two particles are at their 
closest approach in the present time step, causing the algo- 
rithm to scale as O (A/^) . These methods have been reformu- 
lated in terms of DSMC by Lang et al. ifsoll . We typically use 
A/t = 10^ — 10^ test particles, and, by the various improve- 
ments we describe below, in most cases considered, here, we 
can obtain accuracy to within 1%. 

As pointed out in Sec. Ill Al the Boltzmann equation has a 
simple interpretation in terms of the flow of atoms through 
phase space. Hence, the collisionless evolution of the test par- 
ticles is performed by solving Newton's laws for the potential 
U (r, t), and collisions are governed by the collision integral 
[Eq. (O]. The collisions are implemented probabilistically 
(see Sec. IIII C 3b using a scheme that requires the particles 
to be binned into a grid of cells in position space. This serves 
two purposes: (i) It allows for the sampling of the distribution 
function, and (ii) it establishes a computationally convenient 
mechanism for determining which particles are in close prox- 
imity. Thus, the accuracy of DSMC depends on the discretiza- 
tion of the problem, the cell size, the time step, and A/t- It has 
been shown to converge to the exact solution of the classical 
Boltzmann equation in the limit of infinite test particles, van- 
ishing cell size, and vanishing time step llsill . 

In the original DSMC algorithm llssll . a test particle may 
collide with any other particle within the cell. This coarse 
grains position and momentum correlations, such as vorticity, 
to be the length scale of the cells, as observed by Meiburg 
lf52ll . If the cells are not small enough, this transfer of infor- 
mation across a cell could lead to nonphysical behavior To 
combat this, we have employed a nearest-neighbor collision 
scheme 03] outlined in Sec. IIII C"4l where the collision part- 
ner of a particle must be chosen from the nearest neighbors. 
Although nearest-neighbor collisions alleviate this problem, 
the cell sizes still must be small in comparison to the local 
mean-free path and the length scale over which the density 
varies for accurate sampling. 

The time step of the simulation must also be small in com- 
parison to the smallest local mean-collision time to ensure the 
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validity of the basic assumption of DSMC and that physical 
atoms do not propagate further than the local mean free path 
before colliding. To ensure this (and for added efficiency), we 
implement locally adaptive time steps ifssll where, instead of a 
single global time step, the time step can vary over the whole 
system, adapting to the local environment. 



C. Implementation of DSMC 

Here, we consider the basic implementation of DSMC; a 
collisionless evolution followed by a collision step where test 
particles are binned in position space and collisions between 
them are implemented stochastically via a collision probabil- 
ity. We also discuss the various adaptive schemes we employ 
for better accuracy and efficiency, while retaining the desired 
linear scaling of the computational complexity with test parti- 
cle number 



1. Collisionless evolution 

The collisionless evolution is performed by a second-order 
symplectic integrator ifisllsill . which updates the phase-space 
variables of the i*'' test particle in three steps: 



(a) Ax 



{t) 



At 
2rn 



(7a) 



p, {t + At) = p, (t) - At Vq. U{q,,t), (7b) 
At 



r; (t + At) = q, 



2m 



Pi {t + At) 



(7c) 



Symplectic integrators have the properties of conserving en- 
ergy and phase-space volume over long periods of time. 



2. Master grid and LACs 

To perform collisions, we must first bin the test particles 
into a grid of cells according to their position. Collision part- 
ners are then selected from within each cell. In general, the 
binning occurs in up to two levels: (i) the master grid on 
which each master cell is a rectangular cuboid of equal size 
[see Fig. |2ja)] and (ii) the adaptive subdivision of the master 
cells into smaller cells dependent on the number of particles 
in the parent master cell, i.e., LACs [see Fig. |2b)], which is 
an optional refinement. The use of several LAC schemes in 
DSMC is discussed in Ref. f3^ . It is a useful refinement to 
the algorithm for applications to cold-atom systems because 
these typically have large variations in density (such a scheme 
has been employed in Ref. ifssll to account for the large change 
in density during evaporative cooling of a cloud of cesium 
atoms). We now discuss these levels in further detail. 

At the beginning of the collision step, the grid of master 
cells is chosen to ensure all particles are held within its bound- 
aries [see Fig.|2a)]. We choose to keep the size of the master 
cells in each direction constant in time so that if the particles 
spread out further in space during the simulation, we add extra 
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FIG. 2. (Color online) A two-dimensional schematic of the cells used 
for a swarm of test particles, (a) The rectangular master cells are 
all of the same size and are chosen to ensure all particles lie within 
the boundaries of this grid. Cell boundaries are indicated by lines, 
and particles are indicated by dots, (b) An enlargement of two mas- 
ter cells showing their adaptive subdivisions into smaller cells. The 
number of subcells is determined by the number of particles within 
the master cell. 



cells rather than changing the size of the cells. The particles 
are then binned into these master cells, and the number of par- 
ticles in each cell Nc is stored. 

For adaptive subdivision, each master cell is considered in 
turn, and the particles are binned further into a grid of smaller 
subcells according to Nc [see Fig.|2jb)]. Because the number 
of collisions within a cell increases with density (i.e., number 
of particles), the subdivision of highly occupied master cells 
gives a finer resolution of spatial regions where the local col- 
lision rate is highest and, hence, more accurate simulations. 

Our subdivision procedure aims to produce cells in which 
the average number of particles is close to some threshold 
value Nth for which the choice of is discussed in Sec. IIV B 2! 
In our algorithm, we do this by finding the integer / such that 
Nc/2^ is closest to, but not less than. Nth- The master cell is 
then subdivided into 2' subcells, while giving no preference 
to any direction in this subdivision. We choose this division 
scheme over more complicated schemes, as when additionally 
implementing LATS, the protocol for dynamically changing 
grids becomes simpler 

We have adopted the notation of specifying quantities per- 
taining to a particular cell by a subscript c. In what follows, 
when referring to cells, we will mean finest level of cells, i.e., 
the subcells (if used) or master cells otherwise (e.g., if LAC 
is being used, Nc refers to that number of atoms in the sub- 
cell). We do not explicitly label the cells, indeed, this is to 
partly emphasize that the calculations performed in each cell 
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are independent of other cells. Thus, the algorithm is intrin- 
sically parallel and is suitable for implementation on parallel 
platforms (e.g., see Ref. ll4lll ). 



3. Collisions: Scaling 



4. Collisions: Nearest-neighbor selection of partners 

We employ a nearest neighbor collision scheme to com- 
bat discretization effects from finite cell sizes, in particular, 
the so-called transient adaptive subcell (TASC) scheme ifssll . 
Simple sorting of the test particles for the nearest neighbors 



The collision probability for a pair of test particles i and j 
in a cell of volume AV"c is given by 



p 



At 



(8) 



where a (vr) is the total cross section. This collision probabil- 
ity can be derived from the collision integral (l3) via the Monte 
Carlo integration ifisl IstIi . the kinetic arguments [35], or the 
elementary scattering theory ||49| . The correct collision rate is 
established by testing Mc = Nc {Nc — 1) /2 collisions in the 
cell (see Appendix lAl for justification of this choice). This is 
inefficient as the number of operations scales as J\f^, and the 
collision probability may be far less than 1 . However, within 
a cell, the collision probabilities and the number of tested col- 
lisions can be rescaled by a single parameter A such that the 
number of operations scales as A/r llssll . 



P,.. 



A ' 



(9a) 
(9b) 



and still converge to the same Boltzmann equation evolution. 
Here, A is chosen to be 



A 



Mca-^ [Vra (Vr)]^ 



Mr 



(10) 



where [vrtr (wr)],nax '■^^ maximum of this quantity over 
all pairs of particles in the cell and [x] denotes the ceil- 
ing function. This corresponds to Bird's proposal of using 
A = ma.x {Pij} 1I35I1 . while we ensure that Mc is an inte- 
ger and at least one collision is tested (Fig. [T4l demonstrates 
the reduction in collisions if this is not taken into account). 
With this choice of scaling, the maximum collision probabil- 
ity within the cell is < 1 (expected to be close to 1), and the 
number of collisions that need to be tested is reduced to 



Mr 



Nr 



1 



where 



nc ^ aNc/AVc, 



(11) 



(12) 



is the density in the cell. 

This enhancement of efficiency is often missed by other 
stochastic particle methods, or the collisions are adjusted in 
some other manner For example, Tosi et al. [27] introduced a 
scheme for fermions where collision pairs with small classical 
collision probability were neglected. 
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FIG. 3. (Color online) A two-dimensional schematic of how colli- 
sions are perfoimed within the TASC scheme. A single cell (outer 
boundary line) and the distribution of test particles (black dots) are 
shown in (a) and (b) for two different random collisions. The finer 
grid of internal lines represents the boundaries of the TASC subcells. 
The first particle of the collision pair is selected at random from all 
the particles in the cell. In (a), the first particle occupies a TASC 
subcell that contains other particles, and the second participant in the 
collision is chosen at random from these other particles. In (b), the 
first particle (which occupies the central subcell) is the sole occupant 
of a TASC subcell. In this case, we check to see if there are any par- 
ticles in layer I, and if so, the collision partner is chosen at random 
from these other particles. If there were no particles in layer I, we 
would then check layer 2, and so on. 



scales quadratically with the particle number The TASC sort- 
ing scheme retains linear scaling, but it does not guarantee the 
exact nearest neighbor 

The basic TASC scheme is to further bin the particles into 
subcells within the cell [see Figs. [3a) and|3lb)], the number 
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of which is roughly equal to Nc- In our case, the number 
of subcells in each direction is equal and is given by [v^lVcJ 
(with [x\ as the floor function). When a particle is randomly 
picked for a collision, its collision partner is established by 
looking within the immediate TASC subcell [Fig. [3 a)], and if 
not found [Fig.[3lb)], each layer starting closest to the particle 
is searched for other particles. If a layer contains more than 
one particle, the collision partner is randomly chosen from 
that set to avoid any biasing. This reduces the distance be- 
tween colliding pairs significantly and may be decreased even 
more by increasing Afx- 

We use this procedure to select each of the Mc pairs of par- 
ticles for testing if a collision occurs. We also ensure a particle 
does not undergo a second collision in the same time step. 



J. Collisions: Testing and implementation of collisions 

For each of the pairs, the collision goes ahead if i? < , 
where i? is a random number uniformly distributed between 
and 1. As Eq. ([T]i describes binary collisions of point like 
particles that conserve total energy and momentum, only the 
momenta are changed by keeping the total momentum con- 
stant, and the relative momentum vector is rotated about its 
center js^. The scattering angles are determined by using an 
acceptance-rejection Monte Carlo algorithm for the differen- 
tial cross section. 



6. LATSs 
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FIG. 4. (Color online) An example of the sequence of steps in a 
DSMC evolution, (a) A simple DSMC scheme where the whole sys- 
tem evolves according to a single global time step At. (b) An exam- 
ple of a cell using LATS. In this example, the global time step (5t) is 
held constant, while the local time step (Stc) is shown to vary. Col- 
lisionless evolution occurs at each global time step. A collision step 
is performed at the global step when, at least, Stc has passed since 
the last collision step. At global time t^, we show a collision step, at 
which the local time counter (tc) is updated and a new local time step 
(St'c) is established. Here, the local time step decreases, showing two 
further collision steps that follow shortly after the first. 



All of the preceding aspects of our implementation of 
DSMC can be performed with the single global time step At 
for all cells such that the evolution of the system is simulated 
at the times tk = k At, with k as an integer At each of these 
steps, the collisionless evolution is performed, then, is fol- 
lowed by the collision step [see Fig.UJa)]. However, if there 
is large variation in the properties over the system, the use of a 
single time step can be inefficient, as it may be much smaller 
than required for low-velocity or low-density regions. This 
has been addressed by a recent improvement to the DSMC 
algorithm JsH), where a local time step was introduced for 
the collision step. Performing the collision step is computa- 
tionally expensive, so this improvement can lead to a great 
increase in the efficiency of calculations. 

With the use of LATS, there are two time steps of impor- 
tance for each cell: (i) The global time step Sti, which is 
the fundamental increment of time in all cells of the system. 
The global time after k steps is specified as tg = X]i=i ^^i^ 
and during each increment of Sti, collisionless evolution is 
performed [i.e., Eqs. (fTat— (TTcli with At 6ti]. (ii) The 
local time step for the cell Stc, which is the desirable time 
scale for performing collisions in this particular cell. Note 
Sti = nim{Stc}, i.e., we choose the global time step to be the 
smallest value of Stc over all cells in the system at the end of 



each step0 

A collision step is performed at the global time step when, 
at least, a time of Stc has passed since the last collision step 
for the cell under consideration [see Fig. Sb)]. To implement 
this, we introduce a cell timer tc, indicating the time up to 
which collisions have been accounted for in the cell. In gen- 
eral, tc < tg and is incremented by Stc during each colli- 
sion step. Performing collisions in this way ensures that tc is 
within Stc of tg at all times0 and at the end of the simulation, 
all tc are updated to the final time by performing collisions 
with Stc = tg — tc- 

In our simulations, Stc is chosen to be small compared to 
the relevant collision and transit times of the cell. In detail. 



^ If Sti is sufficiently large that the accuracy of the collisionless evolution is 
compromised, Sti is split into smaller increments for this evolution. 

' If a cell becomes unpopulated (Nc = 0), tc may not have been updated 
such that tc = tg before the test particles leave the cell, which decreases 
the collision rate. However, Stc is chosen such that this effect is negligible. 
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these time scales, 



IV. TESTS AND OPTIMAL PARAMETERS 



„coll 



{ric [WrCr(Wr)]max}~' ' 

Axc Aye Azc 



(13a) 
(13b) 
(13c) 



are the mean-collision time, the maximum collision time, and 
the mean-transit times of the cell, respectively. These expres- 
sions are evaluated at the end of each collision step, and the 
average speeds {vx,Vy,Vz) are given by averaging over all 
the test particles within the cell, while Vr(J (vr) is the average 
of Vr(7 (vr) ovcr the particles tested for collisions. The cell 
widths {Axc, Aye, Azc) correspond to the cell under consid- 
eration [e.g., Axc is the adaptive subcell x width if the LAC 
is used, and the master bin width (Ax) otherwise]. 
In terms of these time scales, we take 



Str. = 



{'^collT, 



71 max tr 1 

J '/max'' c I "tr ' c / i 



(14) 



where r^coii, Vmax, and ?7tr are constants less than unity. At the 
end of each collision step, Stc is reset by Eq. (fl4l i. Whenever 
6tc is established without performing a collision step, i.e., be- 
ginning of the simulation or when the subcells are collapsed 
or expanded, we take it to be 6tc = min {ijnmxT^^^ , VtrT^^ }- 

For the accurate simulation of dynamics, it is required that 
Stc ^ Tc°^^ as well as Stc ^ t*'- We also require that it 
is unlikely for a particle to undergo multiple collisions in a 
collision step (accounted for by t™^^)- These requirements 
are ensured by the constants rjcow, f]ma,x, and r/ti , which are 
optimized for the desired accuracy. 

Care has to be taken when the LATS is implemented in con- 
junction with the LAC scheme, as the cells can change dynam- 
ically during the evolution (cells can be resized, can be added 
or can be removed). Our procedure for dealing with dynam- 
ically changing subcells is as follows: As each master cell is 
considered in turn, if the number of LAC subcells changes, a 
new layout of LAC subcells must be established. If the num- 
ber of these subcells increases, then each of these new cells 
inherits the tc of the original cell. Alternatively, if the number 
of subcells decreases, then the new cells are formed by merg- 
ing old cells. In general, the values of tc for each of the cells 
to be merged are different, and we take the new value of tc 
to be the largest of these. This requires tc of the old cells to 
be updated to the new tc, thus, collision steps are performed 
within the old cells before merging, using the time difference. 

When the LAC scheme is implemented with small thresh- 
old numbers (e.g.. Nth < 5) and the number of test particles is 
large (A/r > 10^), it can become inefficient to implement the 
LATS in conjunction with the LAC subcells. In such regimes, 
very dense grids of LAC subcells typically arise, for which 
the computational intensity of the LATS and memory require- 
ments become too great. Furthermore, small cell sizes lead 
to excessively small time steps (e.g., r*'' is proportional to the 
cell size), which further reduces the algorithm efficiency. In 
these cases, it is more efficient to implement the LATS for 
the master cells (i.e., only the master cells have a time counter 
and desired time step) and implement collisions in all the LAC 
subcells using that same desired time step. 



In this section, we develop tests relevant to ultra cold sys- 
tems that we use to validate and to explore how to optimize 
the performance of the DSMC algorithm by quantifying the 
effects of the adaptive enhancements. Primarily, we are inter- 
ested in the quality of the representation of the phase-space 
distribution, since this is of fundamental importance for accu- 
rate Boltzmann evolution. In particular, we address the effects 
of increasing the number of test particles and refining the grid 
on collision rates as compared to exact results. 



A. Analytic results 

We develop benchmark analytic results to calibrate the al- 
gorithm against. To do this we consider the equilibrium 
(Maxwell-Boltzmann) distribution function for a harmoni- 
cally trapped gas 



P 

2m 



U{r) 



feq{p,r)=J\fp{l3nuj) cxp|-/3 
where 



(15) 



(16) 



is a harmonic trapping potential and uj — (ojxUjyUJz) ^ . 
The total collision rate is given by 

Here, we have taken the differential cross section to be veloc- 
ity independent to give a total cross section of <jq. Evaluating 
this expression for the equilibrium cloud, Eq. (fTTt . we obtain 



(18) 



As we are concerned with simulating the collisions of equi- 
librium clouds, it will be useful to consider the instantaneous 
distribution, 

fcoii (p, r) = feq (p + Poz, r) + feq (p - pqz, r) , (19) 

which corresponds to two spatially overlapping clouds with 
equilibrium shapes that are traveling with opposite momenta 
±Po along the z direction. The total collision rate for this case 
is 



7^2 

Rcoii ^^^mPuj^ao 



exp [ -Po- 
rn 



(20) 



For small pq, the term in the square brackets scales as 4 + 
^/3pl/m + {pq), showing that, forpo — 0, Eq. ( l20t reduces 
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to Eq. ([T8]l with Mp 2Afp 
Pq, it scales as 2 + y/pn/rjipo 



as expected. 



While for large 
The first term 



corresponds to the intra cloud collisions, while the linear term 
is that of which is obtained for momentum distributions of 
vanishing width, i.e., Dirac 6 functions S {p ± Pqz). 



B. Grid parameters and test-particle number 

To investigate the accuracy with which collisions are 
treated, we compare the numerical collision rate to the ex- 
act values in Eqs. (fTsT l and ( l20l i. To do this, we calculate the 
relative error of the numerical collision rate and examine its 
dependence on the number of test particles and grid refine- 
ment Q 



1. Numerical collision rate 

For the purpose of comparison, we need to extract a colli- 
sion rate from the DSMC representation of / (p, r, t). To do 
this, we evaluate the mean number of collisions in each cell 
over some time 5tc- Hence, in each cell, the mean collision 
rate is 



Rc 



Ma 

iv) 



AVrSU 



(21) 



where (ij) indicates the indices of the Mc selected collision 
pairs in the cell. Thus, the total collision rate for the system is 



(22) 



oils 



By calculating the collision rate in this way, we are, in ef- 
fect, directly performing a Monte Carlo integration for the in- 
tegral ( fTTl i. which is the basis of the derivation of the collision 
probability in Refs. iflsiljTll . The time step for the cell Stc is 
somewhat arbitrary, and we choose it to give Mc = [Nc/2\ 
collision pairs. 

A convenient length scale for the trapped system is given 
by the thermal widths Wx = ^J^kj^TjmLjl., etc., and we 
choose the master cell widths such that the resolution in each 
direction (relative to these widths) are the same, i.e.. 



Ax 



Ay 



Az 



(23) 



In what follows, 7 will serve as an important parameter to 
specify the fineness of the spatial resolution. 



2. Accuracy 

To increase the accuracy of our numerical calculation of the 
total collision rate, we must improve the accuracy of our rep- 
resentation of continuous distribution / (p, r, t) or take more 
samples. In DSMC, / (p, r, i) is represented in two ways: 
(i) the test-particle swarm, (ii) the grid of cells that sample 
the test-particle swarm. Appendix lAl shows that, without cell 
adaption [i.e., LAC or LATS], (i) and (ii) are largely decou- 
pled. However, simply decreasing the size of the master cells 
can cause large statistical fluctuations in the number of colli- 
sions, as single occupation of a cell becomes more common, 
hence, requiring a larger number of samples. 
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FIG. 5. The relative error of the total collision rate, Eq. (I18t . for 
the equilibrium distribution /eq (p, r) against 7 with A/r = 10^ is 
shown for the cases without (solid line) and with cell adaption where 
Nth = 2 (dotted line), 150 (dashed-dotted line), and 500 (dashed 
line). The results shown here are averaged over 200 initial condi- 
tions, while the error bars give the standard deviation. Without adap- 
tion, the error increases with increasing 7, since /eg (p, r) becomes 
more coarsely grained. However, with the inclusion of adaption, this 
behavior is combatted as the LAC subcells adapt accordingly. We ob- 
tain the initial conditions for the test particles from /eg (p, r) using 
the Monte Carlo acceptance-rejection method. System parameters: 
The harmonic potential is chosen to be the same as that used in the 
ultra cold collider experiment with i^a; = tjy = 27r x 155 Hz and 
cj^ = 27r X 12 Hz, AAp = 2 X 10^ and T = 600 nK. 



Our LAC scheme essentially establishes a local maximum 
size of the cells (i.e., maximum error), which is set by A/r, 
Nth, and n (r, t). In our results, this is seen for the collision 
rate of the equilibrium cloud given in Fig. |5] These results 
show that the magnitude of the relative error does not con- 
tinue to increase with increasing 7 (as it does in the unadapted 
case) but tends to a constant dependent on Nth- With decreas- 
ing Nth, smaller cell sizes are achieved, hence, lower error.0 



The relative error in the collision rate is independent of Mp and ctq . 



^ Care needs to be taken with other adaptive schemes, since the approach 
outlined in Appendix Ia] to remove statistical biasing, neglects to take into 
account statistical fluctuations from other sources (e.g. volume), which 
may become important BtIi . 
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FIG. 6. (Color online) Schematic of the ultra cold collider used in 
Sec. lIV Cl Two clouds initially separated by a distance of 2ro collide 
at a relative momentum of 2po . The number of atoms that have scat- 
tered out of the clouds, after they have passed through each other, is 
referred to as TVs c- 



However, we restrict ourselves to Nth > 2 to avoid the in- 
creasingly large statistical fluctuations mentioned earlier The 
results in Fig. |5]remain qualitatively similar for different val- 
ues of Mt, however, the fluctuations (i.e., error bars Fig. |5) 
increase with decreasing test-particle number 

It is worth noting that systems with identical density distri- 
butions are coarse grained in the same fashion (provided A/r 
is the same when using the LAC scheme), hence, they have the 
same accuracy. For example, the equilibrium (fT5t and colli- 
sion ( fT9l ) distributions have identical relative error profiles as 
seen in Fig. |5] However, if a system is dynamically chang- 
ing and no adaption was employed, evolving to a more dilute 
system would decrease the magnitude of the relative error, 
while increasing if becoming denser For adaptive schemes, 
this is not an issue, as the cell sizes automatically adjust to 
this change. 



3. Performance considerations 

The results in Fig.|5]show that the following cases approx- 
imately have the same relative error in collision rate: [SIMl] 
an unadaptive simulation with 7 = 0.02, [SIM2] a LAC sim- 
ulation with Nth ~ 2 and 7 = 0.2 (we also include the LATS 
for dynamics in SIM2). 

A fuller picture of the merits of using either of these ap- 
proaches for a simulation requires us to understand their re- 
source requirements. 

Speed: We find that, with our code SIM2 is approximately 
five times faster than SIMl for near-equilibrium evolution. 
Note, we only use the LATS scheme in SIM2 for the mas- 
ter cells (as discussed at the end of Sec. IIIIC6I ). It should 
also be noted that this performance indicator is dependent on 
the code implementation and physical problem under consid- 
eration (i.e., equilibrium cloud versus highly nonequilibrium 
situation). 

Storage: SIMl requires ~ 5 x lO*" master cefls, while SIM2 
requires ^ 5 x 10** master cells with a maximum of 4096 
LAC subcells within a master cell (typically requiring a total 
of ~ 7 X 10^ LAC subcells). 



C. Collisions between clouds: Comparison to simple methods 

In this subsection, we consider the collision of two equi- 
librium clouds in a harmonic trap, /e^ (p ± po^, r =F roz), 
shown schematically in Fig.|6] We study this collision using 
our DSMC algorithm and compare its results to a simplified 
model that has been used previously to analyze this problem. 
Initially, the two clouds are centered at locations separated 
by a distance of 2ro along the z direction, chosen to ensure 
that (initially) the clouds do not overlap. The clouds approach 
each other, moving at a relative momentum of 2po, and when 
they overlap, collisions scatter atoms out of the clouds. Here, 
our main interest is the total number of such scattered atoms 
Afsc, after the two clouds have completed passing through 
each other 

The sim ple model we consider was used in Ref. ll32ll (see 
also, Ref. 115811 ) and was derived from the Boltzmann equa- 
tion description of the colliding clouds by making the follow- 
ing approximations: (al) the harmonic potential is ignored 
(collision taken to be in free space); (a2) the momentum dis- 
tribution of each cloud is replaced by S (p ±poz); (a3) the 
dynamics of scattered atoms are neglected. These approxima- 
tions lead to an equation for the densities rii of cloud i = 1 , 2 
of 



d Vr d 

dt '2'd'z 



Hi (r, t) = -Vraani (r, t) n2 (r, t) , (24) 



where Vr ~ 2pi^/m. We can solve these equations directly 
using a pseudo-spectral method. 

An analytic expression may be derived with an additional 
approximation: (a4) The loss of atoms is small enough such 
that the shape of the densities do not deform but remain Gaus- 
sian while the normalization of each cloud A/p decreases. Us- 
ing this, one can integrate Eq. (l24l i over all position space to 
find the total number of scattered atoms from the collision. 



7^2 

Use = —^mPLOxLUyUQ. 

47r 



(25) 



Following the terminology established in experiments, we 
characterize the collider kinetic energy in temperature units 
by the parameter Tcoii = fJ'v'^/2kB, where ji — to/2 is the 
reduced mass. As shown in Sec. IIV Al when considering the 
limiting behavior of Eq. (l20l i. the approximation (a2) is satis- 
fied when Tcoii ^ T (which is the case for collider velocities 
we consider here). That the momentum distributions can be 
replaced with Dirac 6 functions is consistent with many-body 
quantum statistics not playing a significant role in the scat- 
tering that occurs when the two clouds collide. However, the 
internal motion of each cloud can be influenced by quantum 
statistics. 

As the full DSMC solution includes the dynamics of scat- 
tered atoms, it is useful to split the scattered atoms into two 
groups: (i) scattered atoms that have not undergone any subse- 
quent collisions, (ii) scattered atoms that have undergone ad- 
ditional collisions, including all collision partners]^ All of the 



' We include atoms that are scattered out of cloud 1 or 2 by a collision with 
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scattered atoms predicted by Eqs. (l24l and (l25t are of group 
(i). 

In Eqs. (l24l) and dZSl l. A(sc is independent of the details of 
the differential cross section (only depending on the total cross 
section), and this is largely true for the full solution in the case 
considered here. Thus, it is convenient to take (Tq = 87rasc, 
which is of the form of the total cross section for ^-wave scat- 
tering in the low-collision energy limit with scattering length 
ttsc- Additionally, Afgc in both equations is independent of Vr, 
i.e., TcoH- However, this is not the case for the full solution, 
since the collision occurs in a trap. For example, if the ra- 
dial confinement is tight, then a scattered atom can oscillate 
out and back in the radial plane and can recollide (depending 
on the timescale over which the collision proceeds). Here, we 
choose to operate in a regime where these effects are small and 
the simple model should accurately describe the full solution. 
To do this, we choose parameters such that Tcoii ~ 300 /xKQ 
giving a short time scale for the collision. 
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FIG. 7. (Color online) The fraction of scattered atoms due to the 
collision of two equilibrium clouds as a function of a^c. Equation 
i25l [solid green (gray) line] has poor agreement with the solution of 
Eq. ( I24t (solid black line) for TV^c/A/p > 0.05, since approximation 
(a4) is no longer valid. Group (i) scattered atoms (dashed black line) 
and total scattered atoms [groups (i) and (ii)] (dashed-dotted black 
line) from the DSMC solution. The system parameters are given in 
Fig. |5l and for the DSMC simulation, 7 = 0.2, JVt = 10^, and 
Nth = 2. The standard deviation error is not shown as it is on the 
order of the line width. 



The results of Eqs. (l24l i and dZST l, as well as the full solution, 
are shown in Fig.|7]for varying age- All models agree well in 
the low scattering regime Msd-N'p < 0.05, while for higher 
scattering fractions, approximation (a4) becomes invalid, and 
the dynamics of the scattered atoms becomes increasingly im- 
portant. However, the solution of Eq. (l24t agrees to within 
10% of the relative error to the total [groups (i) and (ii) com- 



bined] scattered fraction given by the full solution over the 
whole range. 

We can modify the collision problem and DSMC method to 
a regime that is exactly described by the simplified equation 
(l24l) . To do this, all particles are taken to have momentum 
zkpo along the z axis (the components of momenta in the xy 
plane are zero) and evolve without an external trapping poten- 
tial. Consistent with the approximations going into Eq. ( l24b . 
whenever a pair of particles undergoes a collision, it is re- 
moved from the system (eliminating any need for considera- 
tion of multiple collisions). Due to the form of the distribution 
function, nearest-neighbor collisions cannot be used@ 
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FIG. 8. (Color online) The relative error of J\fsc as calculated by the 
DSMC solution of Eq. J24l l [see text] to that of our pseudo spectral 
solution of Eq. ( 124b . Here, we show the two cases = 10^ (black) 
and 10^ [green (gray)] with 7 = 0.2, Nth ~ 25, and the system 
parameters given in Fig.|5] We use rycou ~ 0.01, T^max = 0.1, and 
rjtr = 0.01. The results shown here are averaged over 200 simula- 
tions, while the error bars give the standard deviation. 



The relative error of Afsc as calculated by the DSMC so- 
lution to that of our numerical solution of Eq. (l24l i is shown 
in Fig. [8] for the two cases A/r = 10^ and 10^. The excel- 
lent agreement of the two results is a good test that the DSMC 
method is correctly implemented. The error bars represent 
the statistical fluctuations of the DSMC results. These fluc- 
tuations reduce with increasing a^c as Afsc increases, while 
between the two cases, they are reduced by a factor of 10, 
since they also decrease with increasing A/r (to be definite, 
these fluctuations are given by the inverse square root of the 
number of scattered test particles). 



an already scattered atom. 
^ For the full DSMC solution, the clouds accelerate as they approach the trap 
center, and we take the value of po that they obtain at the trap center as the 
value to compare against the simple model. 



Particles have no transverse momenta, thus particles from the same cloud 
never leave the proximity of each other Hence, it is required that a par- 
ticle from one cloud is closest to a particle from the other cloud before a 
collision can occur, which results in a decreased number of collisions. 
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V. MANY-BODY SIMULATION OF AN ULTRA COLD 
COLLIDER 

In this section, we demonstrate the application of our 
DSMC algorithm to the simulation of the ultra-cold-atom col- 
lider reported in Ref. ll32ll . The main extension, over the 
DSMC collision test presented in Sec. IIV CI is the inclusion 
of the full two-body collisional cross section needed for a re- 
alistic microscopic description of the collisional interactions. 
We then extend our consideration to the long-time dynamics 
of the collider and how the system progresses to equilibrium. 



A. Collisional cross section 




Tcoii (mK) 



Experiments realizing the ultra-cold-atom collider were 
conducted with *^Rb, which is bosonic, prepared in a single 
hyperfine spin state {F = 2,'mp = 2). The wave function for 
two such colliding atoms is required to be symmetric, hence, 
only the even partial-wave terms in Eq. (|4) contribute to the 
differential cross section. At the collision energies of the ex- 
periment, only the first two even terms contribute, / = and 
I ^ 2 (s and d waves). Thus, the differential cross section 
reduces to 



da 



4 sin^ 5q + 25 sin^ 82 (3 cos^ 



+ 20 cos [5q — 62) sin 5q sin 62 (3 cos^ 9 — l\ 



s- and fi^-wave interference 



(26) 



Taking care to integrate over only half the total solid angle to 
avoid double counting, the total cross section cr ) is given 
by the sum of the total s- and (i-wave cross sections. 



327rft2 



(sin^ (So + 5 sin^ 52) 



(27) 



Calculation of the collision energy dependence of the phase 
shifts 5q and 82 is a nontrivial task. The values that we use 
in our simulations [Fig.|9ja)lare those calculated by Thomas 
et al. and reported in Ref. 0211 . Over the range of collision 
energies shown in Fig. |9]the interference between s- and d- 
wave scatterings can be important, and a cZ-wave resonance 
also occurs. The li-wave resonance can be seen in Fig. |9jb) 
by the peak of the total cross section, attributed to the large 
(i-wave cross section. 



B. DSMC simulations 

Using the full energy and angular-dependent-scattering 
cross section, our DSMC method can provide an ab initio 
prediction for the collider experiments. The full and detailed 



FIG. 9. (a) Numerically calculated i-wave (dotted line) and rf-wave 
(dashed line) phase shifts of Ref. fs^]. (b) 5- wave (dotted line), d- 
wave (dashed line), and total (solid line) cross sections. 
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FIG. 10. Column densities at time -k /2u)x after the clouds reach the 
center of the trap. Tcoii = 200 /iK [(a) and (c)] is a regime of s- 
and d-wave interference, while Tcoii ~ 300 /iK [(b) and (d)] is a d- 
wave regime, (c) and (d) only show the scattered atoms. The initial 
conditions for the clouds are chosen as in Sec. IIV CI and the system 
parameters are given in Fig. |5] while the simulation parameters are 
7 = 0.2, JVt = 10^, and Nth ~ 2. The results were averaged over 
200 simulations. Note, we have compared these results to simula- 
tions with Afr = 10^ also averaged over 200 runs, and we find that 
the number of scattered particles and the angular scattering distribu- 
tions agreed to within 1%. 



comparison with experiments and what information this re- 
veals about the two-body collisions are beyond the scope of 
this paper and will be presented elsewhere (although we note 
that the density images shown in Fig. [T] confirm that our ap- 
proach provides a visually good match to the experimental 
results). 

Here, we present the results of column densities calculated 
after two equilibrium clouds f^q (p ± poz, r =F roz) have col- 
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lided for the cases Tcoii = 200 fjK and 300 fiK. Following 
the experimental procedure ll32l . we calculate these column 
densities at a quarter of the radial trap period (7t/2uJx) after 
the clouds reach the center of the trap. At this time the bulk 
of the scattered atoms reach their maximal extent in the ra- 
dial direction. Figures [TOl' a) and [TOl c) show a regime of s- 
and c/-wave interference (Tcou = 200 fiK), while Figs.fTOl'b) 
and [TOl'd) show a cZ-wave regime {TcoU ~ 300 fiK). Clearly, 
the distribution of scattered atoms deviates from the typical 
5-wave halo (e.g., see Ref. (s^ ). 



C. Long-time dynamics: Rethermalization 

The idea of using rethermalization of colliding condensates 
to perform calorimetry has been proposed in Ref. ll60l . how- 
ever, no direct simulations were made of the thermalization 
dynamics. More generally, there has been significant recent 
interest in how a quantum system rethermalizes ll6ll . partic- 
ularly in systems that might be experimentally realized with 
ultra-cold-atomic gases (e.g., see Refs. ll62l l63ll ). To date, 
much of the attention has been focused on integrable or nearly 
integrable systems where numerical solutions are available for 
small samples of atoms (typically JVp < 10^). However, in 
such regimes, thermalization is often inhibited or strongly is 
effected by constraints (e.g., see Ref. ll64ll ) as well as being 
difficult to explore experimentally due to the small atom num- 
ber (or requiring many similarly prepared systems to get a 
good signal). 

Thus, we are motivated to apply the DSMC method to 
model the dynamics of colliding ultra cold clouds well past 
the first collision. As the collisions occur in the trap, the 
clouds will oscillate back and forth, recolliding each time, 
and thus, are provided with the opportunity to rethermal- 
ize. This system is much larger and classical in nature than 
the small quantum systems generically considered for ther- 
malization studies. However, we believe this is an interest- 
ing system: first as a bridge between quantum and classi- 
cal thermalization in ultra cold gases that is practical for ex- 
perimental investigation. Second, such a system might pro- 
vide a unique opportunity to test the Boltzmann equation in a 
regime where the microscopic parameters are precisely know n 
and with well-characterized far-from-equilibrium initial con- 
ditions. Few equations in theoretical physics have evoked as 
much discussion and controversy as the Boltzmann equation 

- particularly in reference to the introduction of irreversibility 

- such a test could be of broad interest and shed light on some 
long-standing issues in statistical mechanics. 

Our first evidence for thermalization comes from examin- 
ing the density profiles of the colliding clouds at times after 
the first collision. Some examples of these density profiles 
are shown in Figs. fTTI al- fTTl d) and reveal that, as time passes, 
the number of atoms participating in the parametric oscillation 
of the mother clouds along the z axis decreases as the colli- 
sions convert the system to a more isotropic form. Indeed, the 
system clearly appears to increase entropy and approaches an 
equilibrium like configuration. 

In order to quantify the approach to equilibrium, it is use- 
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FIG. 11. Column densities illustrating the long time dynamics of 
rethermalization. (a) At time tco:^ = 7.04 after the clouds have 
passed through each other twice. The colliding clouds are still vis- 
ible (density peaks). When the colliding clouds are depleted, the 
system continues to evolve through collective oscillations that are 
illustrated by the images (b) and (c) at tco^ — 18.85 and 19.60, re- 
spectively. The decay of these collective oscillations occurs on a 
slower time scale than the depletion of the colliding clouds, and the 
distribution does not take on the equilibrium distribution until much 
later times as seen in (d) at tcj^ = 500.02. The trap frequencies are 
Lj, = 27r X 50 Hz and ujx = i^y ~ '2ijJz, and each of the initial 
clouds has A/p — 10® and T — 600 nK. We use an isotropic dif- 
ferential cross section with asc ~ 10 nm. The initial separation is 
chosen such that there is insignificant overlap of the clouds. The mo- 
menta are chosen to give Tcoii = 32.4 /iK, giving a final equilibrium 
temperature of T = 6 /iK. Note, for an isotropic trap, the system 
does not completely thermalize without mean-field effects, since the 
breathing mode does not damp fl2ll . 



ful to consider how various moments of the system evolve dy- 
namically. In Fig. [12] we show the envelope of the oscillations 



vl/2 



in the position spread moment 
characterizing the root-mean square of the distance of the par- 
ticles from the trap center [Note the oscillations of this mo- 
ment occur on a much faster timescale and are shown in an 
inset to Fig. [121] These results show that the system rether- 
malizes quite rapidly over the first approximately five trap pe- 
riods. The number of collisions per particle over the first ap- 
proximately three trap periods is shown in the inset to Fig. [12] 
The steps in collision number, which are initially apparent, 
arise from the periodic recolliding of the clouds. However, 
as the system is distributed over modes, these steps smooth 
out. These results show that, during this initial rapid phase 
of rethermalization, atoms experience > 10 collisions, much 
greater than the value of 2.7 often quoted in the literature from 
the study of Wu and foot fill . 

After this rapid thermalization phase, the relaxation to equi- 
librium proceeds more slowly as energy contained within a 
few low-frequency collective modes waits to be damped. We 
find that two modes are dominant on long time scales. Most 
importantly, a mode that oscillates at 2ujx{^ 2ajy) is domi- 
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FIG. 12. (Color online) Envelope of the oscillations (as seen in the 
lower inset) of the root-mean-square of r. The rapid decay of the en- 
velope within the first ten trap cycles is attributed to the depletion of 
the colliding clouds, while the slower decay is the decay of the col- 
lective modes. The upper inset shows the mean number of collisions 
per atom. 
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FIG. 13. Collective oscillation induced by a small contraction of 
the radial trap confinement for a system with = 2-k x 50 Hz, 
ujx = ojy = lOojz, Mp — 10*, and T = 600 nK. The legend gives 
At used, and as this increases, the results for a single run become 
increasingly indistinguishable from the ensemble-averaged result for 
Mt = 10^ 



nated by radial breathing character and is well described (both 
frequency and damping) by the analytic predictions given in 
Ref. lfl2ll . Also, we note that the rate of relaxation is strongly 
dependent on the trapping geometry and collision rate. 

In relation to thermalization dynamics, it is interesting to re- 
visit the role of test particles in the DSMC simulation. In gen- 
eral, increasing Mt has the effect of reducing fluctuations in 
a simulation and, hence, the number of trajectories needed to 
obtain an ensemble average. However, in order to gain a bet- 
ter understanding of typical results (and, hence, fluctuations) 
that might be expected in experiments, it is necessary to take 
Mt = Mp. To illustrate this, we show some results for a small 
amplitude collective-mode oscillation in Fig.[T3]for a system 
with Mp = 10^ and various numbers of test particles. As the 
number of test particles increases, the results become increas- 
ingly indistinguishable from the ensemble-averaged results. 
However, for Mt ~ Mp, the individual trajectory deviates 
significantly. 

We emphasize that our simulations for thermalization in 
this section have been performed for the case of purely s- 
wave scattering. A detailed study of thermalization, includ- 
ing higher-order partial waves (e.g., as the collision energy is 
scanned across the (i-wave resonance), would be needed for 
detailed comparison with experiments in this area but is be- 
yond the scope of this paper. Along these lines, we would like 
to note an interesting interplay between the partial waves that 
has been shown in the study of the thermalization of mixtures 
by Anderlini and Guery-Odelin ll65ll . In that paper, they per- 
formed an analytical study of near-equilibrium thermalization 
of a two-component mixture and showed that the thermaliza- 
tion time (unlike the collision rate) depended on the interfer- 
ence between the scattering partial waves. 



VI. CONCLUSIONS AND OUTLOOK 

In this paper, we have presented a DSMC method for sim- 
ulating the dynamics of nondegenerate ultra cold gases. The 
motivation for our paper was experiments in which two clouds 
were collided at high relative velocity. In order to simulate this 
highly nonequilibrium regime, we have adopted several mod- 
ern enhancements of the DSMC algorithm (i.e., locally adap- 
tive time steps and nearest-neighbor collisions, introduced in 
other fields) but not previously used for cold-atom simulation. 
We have verified that our algorithm is accurate by compari- 
son to a range of analytic results and simplified models. We 
have also provided some benchmarks of the performance of 
our algorithm against traditional DSMC to quantify the com- 
putational efficiency. 

In order to quantitatively describe the collision experi- 
ments, we have included the full energy dependence of the 
s- and d-wave scatterings in the differential cross section. We 
have presented examples of the scattered distributions for the 
regime of experiments revealing the d-wave shape resonance. 
We have also considered the long-time dynamics of the col- 
liding clouds, allowing them to recollide many times in the 
trap, observing how they approach equilibrium. Our paper 
suggests that this might be a fruitful system for future experi- 
mental studies to test the accuracy of the Boltzmann equation 
and to better understand thermalization. 

A future application of this paper will be to produce a com- 
plete dynamical finite temperature theory. Using a simple 
DSMC algorithm, Jackson and co-workers |ll6l - l20ll have al- 
ready implemented the ZNG formalism ll2lll . In the future, 
we intend to perform a similar extension to c-field formal- 
ism ll66ll . Having efficient procedures for evolving the c-field 
equations that describe the low-e nerg y condensed or partially 
condensed part of the system lH^HSfl, the algorithm described 
in this paper provides the basis for an efficient scheme for 
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simulating the high-energy incoherent modes. Another av- 
enue of investigation that we are currently exploring is an ef- 
ficient and accurate way to simulate the quantum Boltzmann 
equation. That is, to include the effects of Bose-stimulated or 
Pauli-blocked collisions. 
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Appendix A: Number of tested collisions 

The Boltzmann equation describes the evolution of the con- 
tinuous distribution / (p, r, t). However, the replacement of 
/ (p, r, t) with a swarm of test particles introduces fluctua- 
tions, which do not correspond to physical fluctuations when 
A/p 7^ Mr- As a result, hydrodynamic quantities are required 
to be obtained from the averages of mechanical variables, not 
the average of their instantaneous values ll69l . 

In these stochastic particle methods, the collisions of test 
particles inherently average the instantaneous values of the 
collision rate. This leads to a biasing of the total collision rate 
when cells have low occupation numbers (see Fig. fT4l i. To a 
good approximation iItoIi . the probabiUty of Nc test particles 
within a cell is given by the Poisson distribution of which the 

variance is equal to the mean, i.e., 5N^ — — Nc = Nc, 
where 6Nc = Nc — Nc- Note that, formally, the correct num- 
ber of collisions to test (given by elementary scattering theory 
and the derivation of the collision probability from the colli- 
sion integral (|3]l via the Monte Carlo integration iflsilSTll ^ is 



^ c 

2 ' 



(Al) 



However, with Poissonian fluctuations in Nc, we get that the 
mean collision rate is i? oc Mc ^ Nc + Nc (but should be 
cx Nc"). Thus, Poissonian fluctuations can become important 
when the number of test particles per cell is small. However, 
the effect of fluctuations from the finite-test particle number 
can be bypassed (e.g., see Ref. liill ') by instead using the num- 
ber of possible pairs of test particles. 



A/5 



Nc {Nc ~ 1) 



which we have employed in this paper 

To understand the difference in detail, we note that the av- 
erage calculated by the DSMC simulation (denoted by the as- 
terisks) for Eq. (lAlt is 



iV2 =Nc +Nc-Pu 



(A3) 



where Pi is the probability of iVc = 1 (as the simulation ig- 
nores cells with Nc — 1, for which no collisions occur, and 
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FIG. 14. Relative error of the numerical total collision rate in the case 
of the equilibrium distribution ( 115) for the choices of Mc, where the 
error bars indicate the standard deviations of 500 averages. Here, 
the original DSMC algorithm is implemented with Nt = 10^ and 
bin parameter ( 123) 7 = 0.2. The choice Af" is seen to diverge for 
low At as Nc in Eq. l lA3b becomes appreciable, which agrees well 
with the theoretically calculated error (Af " theory) using Eq. ([A3}. 
Using Mc removes this divergence, and the error is seen to agree 
well with the expected error for this discretization (Af^ theory). The 
final data set {Mc no fix) demonstrates the error that arises when Mc 
having non integer values after rescaling, is not accounted for (see 
Sec. lIIIC"3] l [i.e., not including the ceiling function in Eq. dl Ib l. The 
system parameters are given in Fig. [5] 



this must be subtracted from the average). While, for expres- 
sion (IA2I) . 



Nc{Nc-l) =Nc 



(A4) 



(A2) 



which gives the correct total collision rate for the physical sys- 
tem as seen in Fig.fT4l 
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